Semiannual Progress Report No. 5 
To 



Aircraft Guidance and Controls Branch/489 
Guidance and Control Division 
National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665-5225 

For the Grant 

PARAMETER IDENTIFICATION FOR NONLINEAR AERODYNAMIC SYSTEMS 

Grant No. NAG- 1-1065 

Period Covered: October 23, 1991 to April 23, 1992 

From 

Allan E. Pearson 
Division of Engineering 
Brown University 
Providence, Rhode Island 02912 

(NASA-CR-190244) PARAMETER IDENTIFICATION 
FOR NONLINEAR AERODYNAMIC SYSTcMS Semiannual 
Progress Report No. 5 t 23 Oct. 1991 - 23 
Apr. 1992 (^rown Univ.) 34 p 


N92-29329 
Unc 1 as 

G3/02 0035759 


Report Prepared By: 

Principal Investigator 
Allan E. Pearson 




Professor of Engineering 
Tele: 401-863-2602 
Date: May 7, 1992 


Executive Officer 
Tele:401 -863-23 19 
Date: May 7, 1992 




Table of Contents 


1. Introduction 1 

2. List of Scientific Collaborators 1 

3. Completed and Continuing Research 1 

3.1 Pan's Ph.D. Dissertation 1 

3.2 Weighted Least Squares and Comparisons 2 

3.3 Modeling the F-18: Phase One Results 10 

4. Future Research 31 

5. Publications and Presentations 31 

List of Figures 

Fig. 1 Performance of the LS Algorithm for Different Noise Levels 6 

Fig. 2 Performance of the WLS Algorithms for Different Noise Levels 7 

Fig. 3 Performance of the AWLS Algorithm for Different Noise Levels 8 

Fig. 4 Performance of the Prediction Error Method for Different Noise Levels 9 

Fig. 5 Block Diagram and the Best Estimated Linear Models 1 1 

Fig. 6 Bode Plots for (j®)> Hqd (j®) “d H ad 0®) ^ 

Fig. 7 Impulse Responses for (s), Hqj (s) and H^Cs) 1 3 

Fig. 8 Comparing the Physical and Model Responses 14 

Fig. 9 Model Sensitivity to a Shifted Data Set 16 

Fig. 10 Comparing H^s) for AWLS, WLS and LS Algorithms 17 

Fig. 1 1 Comparing Hqd(s) for AWLS, WLS and LS Algorithms 1 8 

Fig. 12 Comparing H ft n(s) for AWLS, WLS and LS Algorithms 19 

Fig. 1 3 Performance Tables for Structure Determination of H^s) and Hqd(s) 21 




Fig. 14 Performance Table for Structure Determination of H^fs) 22 

Fig. 15 Modified Performance Tables for the Constrained AWLS Algorithm 23 

Fig. 16 Comparing H^s) Models for the Constrained AWLS Algorithm 24 

Fig. 17 Comparing Second Order H^jCs) Models 25 

Fig. 18 Study of Residuals for the (2,2) H^Cs) Model 26 

Fig. 19 Study of Residuals for the (4,3) H^s) Model 27 

Fig. 20 Study of Residuals for the (2,1) Hq^(s) Model 28 

Fig. 21 Study of Residuals for the (2,0) H^s) Model 29 

Fig. 22 Study of Residuals for the (4,2) H^s) Model 30 





- 1 - 


1. Introduction 

This Progress Report covers the six month period from October 23, 1991 to April 23, 
1992. The completed work during this period is contained in a Ph.D. dissertation of J. Q. Pan 
and two papers. These are summarized below. Work continues on frequency analysis for 
transfer function identification, and a new study has been initiated into a “weighted” least 
squares algorithm within the context of the Fourier modulating function approach. Progress in 
each of these areas is summarized below. In addition, the first phase of applying these tech- 
niques to the F- 18 flight data is nearing completion, and these results will also be summarized 

below. 

2. List of Scientific Collaborators: October 23, 1991 to April 23, 1992: 

J. Q. Pan Graduate Student Research Assistant 

A. A. Pandiscio Raytheon Graduate Student Fellow 

A. E. Pearson* Professor and Principal Investigator 

Y ; shen* Graduate Student Research Assistant 

♦Received partial support under NAG- 1-1065. 

3. Completed and Continuing Research 
3.1. Pan’s Ph.D. Dissertation 

As implied by the title: System Identification, Model Reduction and Deconvolution Filter- 
ing Using Fourier Based Modulating Signals and High Order Statistics, J. Q. Pan s Ph.D. 
dissertation [1] covers a fairly wide range of topics. Some of the results in this thesis have 
already been summarized in previous reports, and we shall only outline the results here. 
Thus, Chapter 2 entitled “Input Persistent Excitation and Model Structure Estimation” was 
discussed in Section 3.2 of Semiannual Progress Report No. 2, and the contents were 
presented in a short paper entitled On Order Determination for Linear Differential Scalar 
Systems” for the 1991 CISS, Johns Hopkins Univ., Baltimore, MD, in March 1991. Parts of 
Chapter 3 entitled “Frequency Analysis Using Short Time Transient Data” have been dis- 
cussed in Sections 3.4, 3.2 and 3.3 of Semiannual Progress Reports No’s. 2, 3 and 4 respec- 
tively, 1 and abbreviated contents were presented in the short paper [2] entitled “Frequency 
Analysis Via the Method of Moment Functionals” for the 1991 IEEE CDC, Brighton, UK, in 
December 1991. 2 Chapter 4 entitled “Schemes for Model Reduction and Parameter 
Identification in the Frequency Domain” is an application to model reduction of our formula- 
tion for frequency analysis (Chapter 3) using the Fourier modulating function approach. Two 

1 Specifically, that part which deals with single input-single output (SISO) systems. We have 
since dealt with the MIMO case (see Section 3.7 of [1]), and this will be included in a full-length pa- 
per to be written in the near future. 

2 Three preprint copies each of the 1991 CISS and 1991 IEEE CDC papers were mailed to P. C. 
Murphy on September 13, 1991, along with three copies of the paper entitled “Explicit Least Squares 
System Parameter Identification for Exact Differential Input/Output Models” which is scheduled to ap- 
pear in the Proc. of the Eighth ICMCM, Pergamon Press, May 1992. 
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sixth order examples are presented using our method, and these are compared with die 
reduced order models obtained via two other techniques. We plan to prepare a paper on thi 
chapter after completing a paper on the frequency analysis results. Chapter 5 entitled Hig 
Resolution Frequency Estimation in the Presence of Noise Using Complex Sinusoidal Modu 
lating Signals” involves the classic problem of estimating the frequencies of supenmpos 
sinusoids corrupted by white Gaussian noise. This study was briefly menuon^ in Secuon lS 
of Semiannual Progress Report No. 4 and was presented as a short paper [3] at the recen 
1992 CISS, Princeton Univ., Princeton, NJ, in March 1992. The presentation evoked several 
inquiries of interest apparently because of its newness to the signal processing audience and 
the favorable simulation results of this method in comparison with the well-known Yule- 
Walker equation approach. A full length version of this paper has been submitted for publica- 
tion to the IEE E Trans, on Signal Processing. Chapter 6 entitled “Deconvolution and Param- 
eter Identification for Noncausal Nonminimum Phase ARMA Systems Using Inverse Cumu- 
lants” deals with noncausal discrete-time ARMA models driven by non-Gaussian random 
inputs possessing nonzero third and fourth order cumulants. 


3.2. Weighted Least Squares and Comparison With Other Techniques 

' We have begun to examine various ways to ameliorate the deleterious effects of noise 
bevond that which is obtained by our deterministic least squares setting of the Founer base 
modulating function technique. We have found that a “weighted” least squares approach 
seems to hold significant promise, at least for linear system models. To describe these resu ts, 
consider the SISO linear system model: 

A(p)y(t) -B(p)u(t) + e{t), 0<t<T 0) 


where (A(p)Ji(p)) are differential operators, i.e., polynomials in p-d/dt, of a priori order n , 
and e{t) represents the effect of modeling errors. Thus far we have considered the weighting 
situation only for the parametric identification problem, i.e., estimating the coefficients in 
(ACp)JB(p )) given the i/o data [u(r),y(0] on [OX], but we shall eventually investigate 
extending the approach to the nonparametric problem as well, i.e., estimating the frequency 
function H (j co)=fi (t co)/A (i to) given transient i/o data on several intervals [t r ,t r +T], 


r- 1,2 ■ ■ N . 

To encompass stochastic models, a first question might be: What is the effect of the 
modulation process on a continuous-time “white noise” process? To answer this, let us refer 
to the complex form of the Fourier modulating functions of order n which were defined by 
Eqs. (2) and (3) in Semiannual Progress Report No. 4, i.e., 

4> m(r) = e'™** (e' 10 * -D n , 0<t<T = 271/coo (2) 

and its equivalent representation (which follows from the binomial expansion) given by 

*=0 

In the above, i=^T, m is any integer which we shall refer to as the ‘modulating frequency 
index’, (&o=2k/T plays the role of a “resolving” frequency, and n corresponds to the order of 
the differential operator model under investigation. Defining a column vector 
6 = colhaj, • • -a n ,b x , ■ • b n l and applying the Modulation Property, which is embodied in 
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Eq. (4) of Semiannual Progress Report No. 4, the modulated version of the model (1) can be 
put into the (complex valued) linear regression equation format: 

7 tf(m ) = Y(m )0 + £(m ), m=0,l,2 • • 

where the row vector of regressors y(m ) is defined by 

y(m) = row[y/(m), • • y%{m )> ' ' Y «( m )l 

and the (y/(m), y/("»)) defined by 3 

y f(m ) = A” (im C0 o ) n ' ; U (m ), y/(m ) = A" (im 0) 0 )"' ; X (m ). 

The residuals for (4) are given by 


(4) 


(5) 


( 6 ) 


e(m) = je(t )<$> m (t)dt. 
o 


(7) 


If we assume that the modeling error function e(t) in (1) is a zero mean Gaussian sta- 
tionary white noise process with covariance Ee(t)e(t+x)=o5(x) then the question posed at the 
start of the previous paragraph has a nice closed-form answer in that the covariance of the 
residuals e(m ) defined in (7) is found to be given by 

(2n)! 


r o 2 


E e(m )e(m +/ ) = 


— (-1)*- 

T («-/)!(«+/)! 


-n<l<n 


( 8 ) 


0 for all / such that | / | >n . 

If the residuals in the time domain were indeed white, then weighting the frequency domain 
residuals e(m) by the inverse of the Toeplitz matrix W deduced from (8), i.e., for all frequen- 
cies m =0, 1-Af, 

f(-l/ £n)! -n<J<n 

J ^ * (n -/)!(«+/)! (9) 

(W }msn+i ~ 1 o for all / such that | / | >n 

would lead to a Gauss-Markov estimate. Although we cannot achieve the minimum variance 
of a Gauss-Markov estimate, due to the statistical dependence of the regressors on the noise 
e(r), weighting by the inverse of W seems reasonable for a first attempt at reducing the vari- 
ance when the measurement noises are white. Employing this weighting and taking into 
account both the real and imaginary parts of the regressors, the cost function takes the form 

j = (y c -r c e)W e - 1 a r c -r c e> m 

where the following notation applies: 

r/ = [ Re (y( 0), • • y(M) ), Im (7(0), • • 7 <M) ) 1 ^ ^ 

= [ Re (Ytf(O), ‘ ) ), Im (Y<J(0), • • y$(M ) ) ] (12) 

with (y ( m), yj(rn) ) defined in (5) and (6), and the composite weighting matrix W c is defined 


3 Refer to Eq. (4), also Fig. 1, of the Semiannual Progress Report No. 4 for further details, includ- 
ing the display of the n ** order finite difference operator A” . 




by 


W„ 


W 0 

0 w 

fc * 


03) 


The one-shot parameter estimate which minimizes (10) is given by 

=(r c 'w , c' , r c r l r e 'H'-'y e . <m) 

The ordinary least squares estimate (the one we have been working with up till now) is 
obtained from (14) by setting W C =I , i.e., 

«<*- (r c T e )''iyy c . (i5> 

Shen has compared these two estimates in Monte Carlo simulations for a second order system 
with white measurement noises (to be summarized below) and indeed the variance of the 
parameter estimates is reduced in using (14) in relation to (15). 

In addition, we are experimenting with an “adaptive” weighting matrix algorithm, which 
uses the inverse of W c as a first guess, in an effort to achieve a ‘generalized’ least squares 
estimate, thereby further reducing the variance. In this case, the estimate is obtained itera- 
tively by solving for a vector- symmetric matrix pair (Q aw i s ,W a wis) a °f vector-matrix 
equations which have the following form: 

=(r e 'H'^r c r 1 r c 'ty^ 1 ,y c (i6> 


^ awls 8 (® awls ) 


(17) 


where the function giQ^is) is the covariance of the residuals as calculated from the model 
relative to the current estimate of the parameters, 0^, while assuming white independent 
measurement noises on both the input and output signals for this calculation. The standard 
relaxation algorithm for solving this nonlinear pair is to first guess W awls =W c in (16) in order 
to obtain a first guess for 0 aH ,i s , substitute this value into the RHS of (17), thus obtaining a 
new value for W^, which is reinserted into the RHS of (16) to obtain the next value for 
Gan-/*, etc. 4 This algorithm works quite well with a suitable stopping rule based upon a 
chosen threshold for change in the estimate of 6^ . Convergence has usually been obtained 
in five to ten iterations for our applications. 


To compare the above least squares formulations, consider the second order system 
y (r >t-3y (r )+8y (t )=5u (r), 0<r<7\ which possesses the transfer function 


H(s) = 


s 2 + 3s +8 


with the output y (r ) corrupted by additive white Gaussian measurement noise. Two hundred 
Monte Carlo runs were made for each of several noise-to-signal ratios. The input was 


4 In keeping with the property from (9) that W and hence W c is bandlimited to n nonzero values 
on either side of the main diagonal, these iterations are carried out with also bandlimited to n 

nonzero values off the main diagonal, where n is the order of the model under consideration. 




-5- 


«(r)=sinf 2 /5 over a 7=105 time interval for each run. The initial conditions were randomized 
for each run, and the [0,7] time interval was discretized into 256 subintervals so that a stan- 
dard 256=2 8 point DFT/FFT could be used in calculating the Fourier series coefficients of the 
i/o data for the least squares problems. 5 * The results are shown in Figs. 1, 2 and 3 for the 
respective estimates, i.e., Fig.(l) for 0^, Fig. 2 for 0 W ^, and Fig. 3 for 0^. Each figure 
shows the estimates for the three parameters against the ideal values: 0i=8, ©2=3 and 0 3 =5, 
as a function of the RMS noise intensity relative to the uncorrupted signal. The mean and 
standard deviations are shown for each parameter estimate based upon the 200 Monte Carlo 
runs at each noise intensity. In terms of variance, there is a roughly 3 fold decrease in com- 
paring the variance estimate based on 0^ verses 0^ , with the variance for being some- 
what between 0 b and 0^ . The bias is also smaller for 0^ than that for either Q wls or 0^ , 
and takes on significant values only at high noise levels. A more dramatic comparison is 
shown in Fig. 4 where the algorithm used was the popular Prediction Error Method (PEM) 
from the Identification Toolbox by L. Ljung in MATLAB. The superiority of the modulating 
function technique (MFT) is unmistakable in comparing Fig. 4 with Figs. 1-3. A major rea- 
son for this is that the MFT does not have to estimate unknown initial conditions. Another 
reason is due to the fact that the MFT is a direct identification technique for continuous-time 
models, while the PEM first estimates the parameters for a discrete-time model then converts 
this to a continuous-time model. Even with the initial conditions fixed at zero and giving the 
PEM the correct initial state for each run, the MFT algorithm still gave better results (figures 
omitted) though the comparisons are less dramatic in this case. Moreover, the MFT required 
less computer time for each set of Monte Carlo runs (about a five to ten fold decrease depend- 
ing on whether the PEM had to estimate the initial conditions or not). 


5 The formula used for these calculations is based on Simpson’s rule and is given by Eq. (19) 

below. 
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Fig. 1 Performance of the LS Algorithm for Different Noise Levels 
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Fig. 2 Performance of the WLS Algorithm for Different Noise Levels 
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PEM with independently random initial state and initial state guess 
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Fig. 4 Performance of the Prediction Error Method for Different Noise Levels 
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3.3. Modeling the F-18: Phase One Results for Longitudinal Flight 

We have begun to model the F-18 aircraft dynamics based upon information communi- 
cated to us from Drs. Vlad Klein and Gene Morelli. Yan Shen is assisting in this project. 
Our initial investigation has focused on the parameter identification problem for the longitudi- 
nal dynamics using SISO linear differential operator models of the form represented in (1). In 
seeking to find the best model within this class, we have to determine the orders of the poly- 
nomials A (p ) and B (p ) as well as their coefficient values. We were given 40.96s of data 
sampled at 50 Hz which means a total of 2250 samples for each signal. Since our algorithms 
use DFT/FFT techniques to calculate Fourier series coefficients, we truncated the data to 2048 
= 2 11 samples in order to utilize a standard length FFT. The 202 extra data points allowed us 
to test the sensitivity of the model to shifts in the data sequence. To be specific about the 
approximations involved, if z(r) is any of the available signals sampled every 8r=.02s on 


[0,7] with samples zj=z(jbt ), j=0,l -N, N =T /5r , then the needed Fourier series 

coefficients Z(m), m=0,l • • M were calculated using Simpson’s rule, the accuracy of which 
is of order (5r ) 4 : 

T 


Jz( t)e 


-imov 


dt = 


2bt 

3 


z 0 +Z N 


N - 1 

+ £ 2 ZjWffJ 

j= 1.3 • • 


N - 2 

£ 

j=2A 


z j w N 


!PJ 


m =0,1 • ■ M (18) 


where w N =e~ i2n/N . Since N is a power of 2, the above sequence of Fourier series 
coefficients is represented in row-vector form by the first M+l components of the standard 
FFT for the sequence inside the brackets on the RHS of (18): 


_ 25r . . 1 o +z n 

Z = —FFT 

3 


, 2zj, z 2 , 2z 3 , z 4 , • • 2 z w _j 


(19) 


An inspection of the various signal spectra reveals that the bandwidth of the models will 
be on the order of 1 Hz or less. Since the resolving frequency is bounded below by the 
length of the available data, i.e., /o=l/T = 024414 Hz or (OqM). 15340 rad/sec, this implies that 
the number M of modulating functions should be on the order of M~3Q or 40, else too much 
noise will be allowed into the least squares regression. This is based on the fact that the 
highest frequency extracted from the data by the modulating functions is Mf qHz, so we 
impose the constraint: Mf q~\ Hz. Notice that this implies a usage of only the lower 2% of 
the available harmonic frequencies in (19) since N =2048. 6 

The longitudinal motion control diagram is shown in Fig. 5 along with the transfer func- 
tions of our best models based on the available signals. The Bode plots and impulse 
responses of the best models are shown in Figs. 6 and 7 for each of the transfer functions. 
The time domain performance of using these models to predict the outputs in comparison with 
the physical signals is shown in Fig. 8. 


6 The implementation of the n‘ h order finite difference operator A", which is needed in calculat- 
ing the regressors (6), actually requires the first M+n harmonics of the Fourier series coefficients of 
the data. Thus, it is the first M+n components of the FFT sequence (19) which is actually retained in 
implementing the algorithm, where n is the order of the model under consideration. 
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F-18 Longitudinal Motion Control Diagram 



models: 

-0.0489s 2 * 0.0244s - 0.0239 zeros: -0.2497+/-0.6533i 

Ha. (S)= 

S 2 + 0.1 335s + 1 .2097 poles: -0.0668+/-1 .0978i 


Hqa (s)= 


-1.31 95s + 0.2457 


s 2 + 0.2465s + 0.0458 


zeros: 0.1862 

poles: -0.1232+/-0.1750i 


-1.1515 

H- (s)= 

S 2 + 0.5834s + 0.1064 


zeros: none 

poles: -0.2916+/-0.1461 i 


Fig. 5 Block Diagram and the Best Estimated Linear Models 
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Bode Plots of Longitudinal Dynamic Models 

Bode plots of model between ETAH and DE 



Fig. 6 Bode Plots for H qd (J to) and 
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estimated impulse response of 



Fig. 7 Impulse Responses for H & ($ ), H qd (s ) and H^is ) 
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Performance of Longitudinal Dynamic Models 

physical and estimated horizontal tail deflection (rad) of F-18 




Fig. 8 Comparing the Physical and Model Responses 
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In addition to the visual comparison, a quantitative measure of how well the models Perform 
is indicated on each graph as the “output signal-to-output error” ratio in decibels, which we 

define by 


S/E = 20 logio . e(t) = y(t)-yit), 0<t<T 


( 20 ) 


where RMS(y) means the root-mean- square value of the output signal y(t) over the [0, ] 
time interval, and y(t) is the estimated output using the model. Since our algorithm does not 
estimate initial conditions, we used a standard Luenberger observer running backwards in time 
(details omitted here) in order to determine an initial condition for each model output y{t) 
when comparing the signals depicted in Fig. 8. 

A test to determine the sensitivity of the above models to shifts in the [0,7*] data sets is 
shown in Fig. 9 for a 200 point-shift. (The shifts in the pole-zero plots are smaller for shorter 
length point-shifts.) The scatter in the pole-zero plots reveal that the H de {s) model is the 
most sensitive of the three models in this regard. 

As implied by the discussion in Section 3.2 above, we have been experimenting with tv\o 
variations on the standard least squares algorithm using the Fourier based Modulating Func- 
tion Technique (MFT). The symbols “LS”, “WLS” and “AWLS” are used to distinguish 
respectively, the least squares estimates given by (15), (14) and (16)-(17). The models 
obtained using these three algorithms are compared in Figs. 10, 11 and 12 for the three su - 
systems. In addition to the visual comparison, the S IE ratios show that the adaptive weighted 
least squares algorithm is the best in each case. 
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Fig. 9 Model Sensitivity to a Shifted Data Set 
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physical and estimated horizontal tail deflection (rad) 



physical and estimated horizontal tail deflection (rad) 



physical and estimated horizontal tail deflection (rad) 



Fig. 10 Comparing ( s ) for AWLS, WLS and LS Algorithms 
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phjsical and estimated pitch rate (rad/s) 





Fig. 11 Comparing H qd (s) for AWLS, WLS and LS Algorithms 
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physical and estimated angle of attack (rad) 



Fig. 12 Comparing H^is) for AWLS, WLS and LS Algorithms 
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Validation of the models with respect to the orders of the numerator and denominator 
polynomials in the transfer functions is the object of the materials presented in the remaining 
Figs. 13-22. The “AWLS” algorithm has been utilized for each of the models obtained in 
these figures. Quantitative data comparing various models for each subsystem is summarized 
in the “performance tables” of Figs. 13 and 14. As indicated in the tables, model orders up 
to and including 5 were tested. In addition to the S/E ratios, the minimum values of the cri- 
terion function J (defined in (10)) is shown in these tables for each pole/zero pair, as well as 
the stability property of the resulting model. Our best transfer functions summarized in Figs. 
5-7 were based upon the data in these tables, though “best” is subject to some qualification. 
For example, using the notation “(2,0)” to imply a model with 2 poles and 0 zeros, and 
“(2,1)” a model with 2 poles and 1 zero, etc., the table for H^is) in Fig. 13 shows that the 
(4,3), (4,4) and (5,4) models are each slightly better than the (2,2) model. However, the 
difference is not deemed significant enough to justify the more complex model and, therefore, 
the (2,2) model is offered as the best for H de (s). A further study of the residuals for the 
H de (s) substem is discussed below. 

An inspection of Figs. 13 and 14 also reveals that unstable models were obtained for 
various structures, sometimes with quite respectable S/E ratios. Based upon the supposition 
that a stable model is preferred, a modification of the AWLS algorithm was tested that forced 
the model for the end result to be stable. 7 The performance tables for this “constrained” 
AWLS algorithm are shown in Figs. 15 and 16. This did not change the conclusion for the 
best models of H de (s) or H qd (s), but the conclusion for the H^s ) is drawn into question by 
comparing the (2,0) and (4,2) models. The Bode and impulse response functions for these 
two models of (s) are shown in Fig. 17 (also, the pole-zero plots). The differences are 
clearly minuscule. However, the predicted outputs for these two H^s) models is compared 
with the physical angle of attack output at the bottom of Fig. 16. Here it is seen that there is 
a noticeable improvement in the (4,2) model over the (2,0) model, but the advantage is prob- 
ably not justified to warrant the additional complexity. Further comments regarding the resi- 
duals for the H^s ) models are given below. 

A study of the residuals for various models of the three subsystems is contained in Figs. 
18-22. The two top graphs in each figure show respectively, the time and frequency plots of 
both the output signal and error signal e =y -y . The bottom two graphs in each figure show 
the modulated residuals (on the LHS), and their normalized covariance (autocorrelation) as a 
function of the modulating frequency index (on the RHS). The X 2 95% confidence limits for 
judging whether residuals are white, or not, are drawn in dotted lines in the lower right hand 
graph for each Fig. 18-22. Based upon this test, each of the models (more or less) can be 
accepted as the “true” model. Again, simplicity weighs towards the lower order model when 
two models satisfy this test as in the case of the H^is ) and (s ) subsystems. 


7 This was accomplished during the iterative solution of (16) and (17) by testing the stability of 
the system for each iterate of 0^ , then starting the iterations over again (if unstable) using as a start- 
ing point the value of 0^ obtained by reflecting the unstable poles of the system about the imaginary 
axis thereby creating a stable system. 
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Performance table of model Hde(s) under different structure assumptions 


zeros 

poles 

0 

1 

2 

3 

4 

5 

i 

•/•■-O.OIdB 

xifeto* 

*>--1«.07dB 

unstable 

>3.41x10* 





2 

«*-0.42dB 

>•87*10* 

*>-31MB 

stable 

X.O0X1O* 

stable 

>722x10* 




3 

t/s— O.lTdB 
unstable 
>•21x10* 

•M-‘107dB 

unstable 

X2SX10 10 

*>«3.56dB 

stable 

>*ja<10’ 10 

*--8.6*0 
unstable 
>2.76x1 O' 11 



4 

70dB 

liable 

>1.62x1 O' 7 

«/#« -8.01 dB 
unstable 
>123x1 0 7 

t/^4.51dB 

stable 

>2.61x10T 7 

*>»1026dB 

■table 

>1.16x10* 

«M«10.16dB 

■table 

>4.42x1CT 10 


5 

s/a*«.20cfi 

unstable 

>12Ex10* 

*/•* -86.97dB 
unstable 
X63X1CT 8 

«/•* -16.40dB 
unstable 
>2.87X10* 

s^«3.07dB 

unstable 

XSftcIO* 

a ^ J tin 

«/*a102SdB 

■table 

>i 2 e<io- 11 

r/AVAJI c 

*»4 03dB 

unstable 

>*.6*<i<r 12 j 


s/e: signal-to-error ratio. J: criterion value, algorithm: unconstrained MFT/AWLS. 


Performance table of model H qd (s) under different structure assumptions 


zeros 

poles 

0 

1 

2 

3 

4 

5 

i 

t/**-2.28dB 

stable 

j-a.afcticr® 

*>■- 226dB 

stable 

X2b(10* 





2 

*>>1.02dB 

stable 

>1.02x1O 10 

1.0MB 
stable 
>2.34x1 CT 11 

t>-1124dB 

■table 

>2.801 O’ 11 




3 

*«1.8MB 

stable 

>4.82x10 10 

«/»-2.82dB 

stable 

X7b<10 10 

*«-89.96dB 

unstable 

M0b«r 12 

*» -88.7MB 
unstable 
X.74xicr 12 



4 

W**1.94dB 

stable 

>2.82x10T 10 

*-0.8MB 

>1.»c10* 

«/*--26.7tkfi 

unstable 

>1.71x10* 

**124dB 

unstable 

X.OtxlO 18 

unstable 

x.aexio ’ 18 


5 

*— 5.36dB 
unstable 
XTIxlO 10 

*--82.51dB 

unstable 

X.17X10- 10 

Wr>-74.18dB 

unstable 

X.06x1O 10 

**-46.17dB 
unstable 
>1.18X10* 
ur 

8.61 dB 
unstable 
J-1.38X10' 11 

T/AYA/l C 

s/#* -17.87dB 
unstable 

xoaxio -11 


s/e: signal-to-error ratio. J: criterion value, algorithm: unconstrained MFT/AWLS. 


Fig. 13 Performance Tables for Structure Determination of 0 ) and H qd (s) 
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Performance table of model H ac j(s) under different structure assumptions 




MdB 

fUbfc 

>2.3^10^ 

i/a*3.04dB 

stable 

>i.*ixicr® 

«/a^2SdB 

>i.66xicr 10 

t/a-0.17dB 

stsblf 

>1.68X1 O' 10 

•/a-4.O0dB 

stable 

J-1 22x1 tr® 

«/a-2£8dB 
unstable 
>7.80X1 or 12 

«/a--123.7dB 

unstable 

j-g-eartcr 11 

i/#~-112.0dB 
unstable 
>3.61x10^ 1 

t/a* -41.31 dB 
unstable 

J-4.81XKT 11 

t/a*-6^4dB 

stable 

>*.01x1 (T 13 



unstable unstable 

>7.87X1CT 12 X.88X1CT 12 


(M-IO.&OdB «fr^. 6 0d B 

unstable unstable 

>i.06xi<r 12 >i.oixier 12 


a/*=-6.06dB 

>7.52x1 (T 13 


s/s«-3772dB 

unstable 

•W.00*1(T 14 


stable 

>7.17*1 cr 13 


«/a«-11.47dB 

unstable 

>1.00*1 or 13 



s/e: signal-to-error ratio. J: criterion value, algorithm: unconstrained MFT/AWLS. 


Fig. 14 Performance Table for Structure Determination of H^is) 
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Performance table of model H(j e (s) under different structure assumptions 


zeros 

poles 

0 

1 

2 

3 

4 

5 

i 

•table 

xiexicr 8 

*-*.36dB 

>*.e 6xi or® 





2 

**-0.42dB 

stable 

>^.87x1 or 8 

•/•-8.1MB 

>*. 08 x 10 ® 

•/•-4 4MB 
•table 

>72&ior® 




3 

-O.OMB 

>*.38X10® 

W#-0.7MB 

•tabic 

>1.1(k10® 

^>*0.66dB 

■tabic 

>2£3X10 10 

stable 

j-a.o&iicr 11 



4 

•M-OTOdB 

•table 

>1.52k10 7 

tM-IJMB 

■tabic 

>1.61x10 7 

«/M4£1dB 

•table 

>2.51x1 or 7 

•/•H02MB 

•table 

>i.ia<io® 

•table 

>4.«a<i(r 10 


5 

■ 


TTlfc/: ■ 

■ 

HiH 

i/#«1023dB 

•table 

>1 .Mxior 11 

W*-9..7ldB 

•table 

>*.96x10 12 


s/e: signal-to-error ratio. J: criterion value, algorithm: constrained MFT/AWLS. 


Performance table of model Hqd(s) under different structure assumptions 


zeros 

poles 

0 

1 

2 

3 

4 

5 

i 

«^*=-2.28dB 

liable 

xjsfcicr 9 

gtmhlc 

^^akior 9 





2 

t/#-102dB 

•table 

>1.02k1(T 10 

t/#*1 1 06dB 
•tiblc 

J^MxKT 11 

*/**ii.24<e 
■table 
>2Mki or 11 




3 

«M-1.82dB 

X92x10 10 

•/•■2.8MB 

>*.74x10 10 

•/•--1.8MB 

xinxio 11 

-1.8MB 

•table 

J4.98XKT 11 



4 

•table 

>2.82x1Cr 10 

•/*-0.83dB 

stable 

>i 9^<icr® 

•/•«2.31dB 

ifthlc 

>2.48x10® 

•/•-4.9MB 

•tabic 

>*33x10 18 

i/»-4.68dB 

•tabic 

>i.o8xior 12 


5 

*/*-2.45dB 

stable 

x.iacio 10 

t/#=-4 56dB 

iUhk 

>*.34xio® 

«/•* -8.7MB 
•tabic 

>4.08x1 cr* 

•/•«-£ 44dB 
>*.78X10- 

•/•-25MB 
>2.37x1 CT 10 

•*•*3.01 dB 

stable 

>2.47X1 or 10 


s/e: signal-to-error ratio. J: criterion value, algorithm, constrained MFT/AWLS. 


Fig. 15 Modified Performance Tables for the Constrained AWLS Algorithm 
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Performance table of model H ac j(s) under different structure assumptions 


zeros 

poles 

0 

1 

2 

3 

4 

5 

i 

*~1.94<S 

•table 

x.«Kicr* 

W>-3.e<dB 

stable 

>1.*1xia* 





2 

•table 

«f**.17UB 

liable 

>i.e8xicr 10 

i/*-454dB 

•tabic 

J^.Kklor 11 




3 

i*-4.06dB 

liable 

j-i-adcr* 

j^wxicr 11 

Ws- 4.8MB 
J-4.00X1CT 11 

*-1.2MB 

liable 

^.sexier 11 

— 



4 

liable 

>234x1 CT 10 

•/•=-3.86dB 

stable 

>3.1*xicr 10 

«/*-11.91dB 

■table 

> 1 . 1 flxicr 12 

t/fc-10.86dB 

stable 

>i 2 tt<icr 12 

s/*^2&dB 

stable 

>7.17X10T 13 


5 

*/*=-1.73dB 

liable 

J=2.06x1tT 10 

i/s* -6_24dB 
stable 

>e.oixior 13 

t/#*-6.06dB 

stable 

>7.62x1 O' 13 

«/**)75dB 

■table 

>7.72x1CT 13 

t/**-0.66dB 

liable 

x.eexicr 13 

■table 

>2.98*1 CT 13 


s/e: signal-to-error ratio. J: criterion value, algorithm: constrained MFT/AWLS. 



Fig. 16 Comparing H^is) Models for the Constrained AWLS Algorithm 
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(2,2) model for DE-ETAH 


DE and its time domain residuals Spectrum of DE and its residuals 



time(s) frequency (Hz) 


Modulated Residuals Normalized Auto— correlations 



Re modulating frequency index 


Fig. 18 Study of Residuals for the (2,2) H^is) Model 
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(4,3) model for DE-ETAH 


DE and its time domain residuals 



Spectrum of DE and its residuals 



frequency (Hz) 


Modulated Residuals Normalized Auto-correlations 



Re 


modulating frequency index 


Fig. 19 Study of Residuals for the (4,3) Model 
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Fig. 20 Study of Residuals for the (2,1) H^(s) Model 
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(2,0) model for ALFAVG-DE 




frequency (Hz) 



Fig. 21 Study of Residuals for the (2,0) H^U) Model 









Fig. 22 Study of Residuals for 
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Spectrum of ALFAVG and its residuals 



0.00 0.25 0.50 0.75 1.00 

frequency (Hz) 


Normalized Auto-correlations 



0 5 10 15 20 25 

modulating frequency index 


(4,2) H^is) Model 
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4. Future Research 

We shall continue our investigations into modeling for the F-18 aircraft dynamics, focus- 
ing on the lateral dynamics as the next phase. We shall also consider obtaining the frequency 
transfer function models directly for this system using the theory put forth in [2], extended for 
MIMO systems. One further point of interest relates to the model reduction problem, which 
is the subject of Chapter 4 in Pan’s thesis [1]. Shen has just rerun the two examples in Sec- 
tion 4.5 of [1] and found that the AWLS (adaptive weighted least squares) algorithm appears 
to give significantly better lower order approximations in the frequency domain than Pan’s 
results which, in turn, were better than the methods in the literature with which Pan compared 
his algorithm. If this turns out to be generally true, then the weighted least squares formula- 
tion will have even broader utility than we anticipated. 
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